The nature and strength of inter-layer binding in graphite 
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We computed the inter-layer bonding properties of graphite using an ab-initio many body theory. 
We carried out variational and diffusion quantum Monte Carlo calculations and found an equilibrium 
inter-layer binding energy in good agreement with most recent experiments. We also analyzed the 
behavior of the total energy as a function of interlayer separation at large distances comparing the 
results with the predictions of the random phase approximation. 
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The excitement generated by the ability to fabricate 
graphene layers and possibly to tune their electronic 
properties has renewed great interest in weak interactions 
in graphitic systems [l[ . The hope of using graphene as a 
component for next generation electronics relies, among 
other things, on a detailed understanding and control of 
how it interacts with its surrounding (e.g. with support- 
ing substrates) @ . 

Unfortunately the nature and strength of binding in 
graphitic materials are poorly understood. For exam- 
ple, large uncertainties are associated with a fundamen- 
tal physical quantity such as the strength of inter-layer 
binding in graphite. In addition, the way the interac- 
tion between planes decays as a function of distance is 
controversial^] [4] , casting doubts on our current under- 
standing of weak binding in carbon based system. Both 
binding strength and power law behavior of interlayer in- 
teractions in graphite are relevant to the comprehension 
of a multitude of materials, including graphite intercala- 
tion compounds, novel nano-electronic components, and 
carbon based systems for hydrogen storage. 

From a theoretical standpoint, unravelling binding in 
graphitic systems is intimately related to understanding 
the role played by dispersion forces, and to acquiring the 
ability to describe these purely quantum mechanical in- 
teractions at a high level of accuracy. Local (LDA) [5J] or 
semi-local Q approximations to Density Functional The- 
ory (DFT) do not correctly describe long range correla- 
tion, due to the local character of the exchange and cor- 
relation potential. Progress 0] has been recently made 
in including dispersive interactions within a DFT for- 
malism, in a self consistent, non empirical manner and 
binding between graphene layer has been predicted. In 
Ref. [7J the authors report an interlayer distance overes- 
timated by 7% with respect to experiment, and a binding- 
energy of 45.5 meV/atom consistent with most extrapo- 
lated measurements of exfoliation energy [l2| . (Note that 
the calculations of Ref. Q are for two isolated graphene 
layers, not for a graphite solid). Semi empirical methods 
have been often used to treat dispersion correla- 
tions where DFT energies are corrected with a contri- 
bution coming from attractive Cq/R 6 potentials between 



pair of nuclei. However the power law behavior used in 
empirical and semi-empirical approaches to describe non 
retarded dispersion forces has been recently questioned 
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At present, no direct measurement of graphite binding 
energy is available. However estimates based on theoret- 
ical models have been reported in the literature, using 
experimental data for exfoliation energies (EE, the en- 
ergy required to remove one graphene plane from the 
surface of a graphite solid}. Three experiments have re- 
ported data for EE [l(| HI 12] j and a common aspect 
to the experimental analyses is the use of simple, fitted 
force fields to model either C-C or C-H interactions. The 
early work of Girifalco and Lad gave a EE value 
of 43(5) meV/atom. Using a Lennard- Jones potential, 
the difference between exfoliation and cleavage energy 
(that is the interaction between two semi-infinite crys- 
tals) was estimated to be 18% but the exact difference 
remains unknown. The work by Benedict et al. 
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trapolated the interaction energy between graphite layers 
(33 meV/atom) from measurements on collapsed nan- 
otubes, using a force field to model the tubes elastic 
properties. More recently Zacharia et al. [lfj| performed 
detailed desorption experiments of aromatic molecules 
from a graphite surface. The graphite EE was derived 
by extrapolating the molecules EE as a function of the 
number of carbon atoms, thus obtaining a value of 52(5) 
meV/atom. This yields an estimate of the cleavage en- 
ergy (62 meV/atom) which is twice as much as that re- 
ported in Ref. fill ]. 

Given the state of experiment and theory in determin- 
ing the binding in graphite, there is a clear need for ac- 
curate calculations, eliminating as much as possible all 
approximations used so far in the literature, and possi- 
bly providing guidance to future experiments. Here we 
report the binding curve of graphite in AB stacking as 
obtained using quantum Monte Carlo (QMC) calcula- 
tions, that is a many-body computational technique (l3| 
capable of accounting for dispersion forces [3] , [3] • We 
obtain an equilibrium inter-layer bond distance in satis- 
factory agreement with experiment ( 3.426(36)A versus 
3.35A ), and a binding energy [29| of 56(5) meV/atom, 



in accord with the mesaurements of Zacharia et al. [l(| ■ 
We find that at distances between 4 A and 8A the total 
energy curve exhibits a ~ D~ 4 2 behavior as a function of 
inter-layer spacing D, which is very similar to that pre- 
dicted by the Random Phase Approximation (RPA) ap- 
plied to two-dimensional (2D) semi-conducting layers[3|- 

In our investigation we have carried out variational 
Monte Carlo (VMC) and Lattice Re gula rized Diffusion 
Monte Carlo (LRDMC) calculations [ljjjj with the Tur- 
boRVB code [17j. Our many body wavefunction is the 
product of a Slater Determinant and a Jastrow many 
body factor. The determinant is obtained with N/2 
molecular orbitals ipj (r) , each doubly occupied by oppo- 
site spin electrons ( N is the total number of electrons). 
The orbitals ipjif) are expanded in a Gaussian single- 
particle basis set {fa}, centered on atomic nuclei, i.e. 

Electron correlation effetcs are included in our 
wave function (WF) through the Jastrow factor 
J(n,-- ■ , r N ) = n exp(/(r 4 , r}), where f(f, f) depends 

i<j 

only upon two-electron coordinates. The function / is 
expanded in a basis of Gaussian atomic orbitals fa,: 
f(f,f r ) = J2i j 9i,j4>i(r)4 > j('r')- The convergence of this 
expansion is improved by adding an homogeneous term 
and a one body contribution, thus satisfying the electron- 
electron and the electron-ion cusp conditions, respec- 
tively [H|][l8j]. The basis set used for the Jastrow in- 
cludes 2s2p Gaussian orbitals. By optimizing the coef- 
ficients c/ij, we can treat in a non perturbative way the 
dynamical transitions to high angular momentum atomic 
states for pairs of electrons localized around each atoms. 
As shown for the benzene 15 1 and water dimer [l9[ , these 
transitions are responsible, at the first order of pertur- 
bation theory, for the weak attractive dispersive forces 
between atoms at large separation. 

In the following, the molecular orbitals ipj in the Slater 
determinant are obtained from a self-consistent DFT- 
LDA calculation. One may then optimize only the Ja- 
trow factor, by keeping fixed the determinant built from 
LDA orbitals (hereafter referred to as J-DFT WF ap- 
proach); alternatively one may simultaneously optimize 
both J and the determinant using the method described 
in Ref. 20]. The minimal Gaussian basis set required to 
build an accurate determinant was chosen by compar- 
ing graphite binding energies (BE) [2i| as obtained using 
plane waves (PW) and Gaussian basis sets (see FigQ}. 
We note that at each atomic positions, PW calculations 
are free of basis set superposition errors and they can 
be converged by controlling one parameter, the kinetic 
energy cutoff. In the case of Gaussian, we used an even 
tempered local basis where the parameters ai and [3i of 
the Gaussian exponent Zi = on(3\ of each angular mo- 
mentum I, were optimized by performing a series of total 
energy LDA calculations. We considered two basis sets: 
4s4p2d and 8s8p4d. For both of them we computed the 
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FIG. 1: Binding Energy (BE) curve for graphite with AB 
stacking obtained at the DFT-LDA level of theory, using plane 
waves (PW) and Gaussian basis sets. The results of PW 
calculations (solid black line) converged as a function of the 
kinetic energy cutoff (90 Ry) are in excellent agreement with 
those of 8s8p4d Gaussian basis sets (red triangles), carried 
out with the same cell (32 atoms and the P point). Smaller 
Gaussian basis sets (4s4p2d, full circle) yield a much larger BE 
and a slightly smaller equilibrium distance. Fully converged 
PW calculations (4-atoms unit cell, 90 Ry, 20x20x8 k-points) 
yield a BE of 24 meV/atom and an equilibrium inter-layer 
position D of 3.30 A . Note that D is significantly affected 
by convergence as a function of k-point sampling, while the 
value of the BE depends weakly on it. 



binding energy (BE) of graphite at the LDA level for a 
system of 32 atoms using only the T point. The 8s8p4d 
basis set reproduces the same BE curve obtained with 
PW converged with respect to the kinetic energy cutoff 
(90Ry). 

In our QMC calculations we simulated a 2 x 2 x 1 and 
2x2x2 super-cell with periodic boundary condition, con- 
taining 32 (128 electrons) and 64 atoms (256 electrons) 
respectively. The carbon valence-core interaction was de- 
scribed by a energy-consistent pseudopontential 2l|. In 
our calculations we fixed the in-plane geometry to the 
one determined experimentally (C-C distance = 1.42 A 
). We verified, at the LDA level, that a change as large 
as 10% in the in-plane carbon-carbon length affects the 
inter-plane equilibrium distance by less than 3%, while 
the BE changes by 1 — 2meV/atom. 

In TableUwe compare the results obtained by full wave 
function optimization with those of the J-DFT WF ap- 
proach for graphite at inter-layer distance D = 3.7 A 
. The latter provides a reasonably accurate variational 
guess. However, full optimization of the WF parameters 
(including the exponent of the basis set) yields a decrease 
of energy per atom of 34(3) meV/atom, which is of a rele- 
vant magnitude to the energy scale we wish to investigate 
in this work. The LRDMC energy is much less sensitive 
to the initial state (guiding function) , used in this ground 
state projection technique. In fact the optimization of 
Jastrow and Determinant in the guiding function, leads 
to a consistent LRDMC energy. The quantitative agree- 



Basis Set 


Method 


E (eV/atom) 


<j(eV/atom) 


4s4p2d 


Opt. VMC 


-154.428(3) 


1.79(5) 


8s8p4d 


J-DFT VMC 


-154.505(1) 


1.66(1) 


8s8p4d 


Opt. VMC 


-154.540(1) 


1.60(1) 


4s4p2d 


Opt. LRDMC 


-154.787(8) 




8s8p4d 


J-DFT LRMDC 


-154.891(5) 




8s8p4d 


Opt. LRDMC 


-154.899(5) 





separation D [ A ] 

5.0 6.0 7.0 8.0 



TABLE I: Variational (VMC) and Diffusion (LRDMC) total 
energy and energy variance (a) obtained using two different 
guiding wave functions for the 2x2x1 super-cell at a sepa- 
ration distance D = 3.7 A. J — DFT denotes a wave function 
with optimized Jastrow and a determinant part from LDA 
calculations. "Opt." refers to the guiding wavefunction where 
both the Jastrow factor and the orbitals were simultaneously 
optimized. 



merit between the VMC and LRDMC calculations con- 
firms that the key ingredients of the electron correlations 
are already included in our Jastrow factor. In the follow- 
ing the results for the 2x2x1 super-cell are obtained 
by fully optimizing the wavefunction, while, due to the 
computational cost of the optimiziation, we will present 
only LRDMC calculations for the 64 atom system. 

After assessing the accuracy of the guiding wave func- 
tion, we considered finite-size (FS) effects. We expect 
the errors due to FS effects in the in-plane directions 
(i.e . x and y directions) to cancel out to a large extent 
14j . as we compute energy differences between systems 
(graphite and graphene) with rather similar bonding and 
electronic properties. In-plane FS errors arising from the 
kinetic and Hartree terms (one-body corrections) can be 
treated within a standard DFT approach with appropri- 
ate k-point sampling. Other FS errors come from the 
artificial periodicity of the exchange-correlation hole due 
to the periodic Coulomb potential. Kwee, Zhang and 
Krakauer (KZK) [23j proposed to calculate the two-body 
corrections within LDA where the exchange and correla- 
tion energy is replaced by the LDA energy parametrized 
for a finite system. We applied KZK corrections as im- 
plemented in the PWSCF code 
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We cannot rely on any error cancellation in the z direc- 
tion, i.e. F.S. errors due to a finite number of graphene 
layers in the simulation cells. In addition, the KZK 
method cannot provide a robust correction scheme in this 
case due to the lack of long range effetcs in the LDA ex- 
change and correlation functionals. We estimate the long 
range behavior of the interaction between planes by fit- 
ting the results of calculations performed on the 2x2x1 
super-cell at distances D > 4 A . These VMC (LRDMC) 
results are reported in Fig. ([2]) and show a behavior 
E{D) ~ D- a with a = 4.2(1) (4.2(3)). Although the 
LRDMC data are affected by larger error bars, we can 
safely conclude that a > 4. We note that using the RPA 
applied to 2D systems, Dobson et al. Q found a power 
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FIG. 2: Binding energy (BE) of graphite as a function of 
separation between planes (a), and of the number of layers 
(n) included in the periodic cell used in our calculations (b). 
Note the logarithmic scales. In (a) results obtained with VMC 
and LRDMC are reported. In (b) only LRDMC results are 
shown. 



law behavior ~ D~ 3 log(D / Do) for infinite 7r-conjugated 
layers. One does not expect to see the asymptotic form 
of Ref. H in our work because the unusual interaction, 
arising from coupling between long-wavelength fluctua- 
tions in the plane, it is expected to arise at much larger 
distances than those studied here [3l[ 

Using the power law determined in our calculation we 
can derive a scaling relation between the graphite BE and 
the number of layers n. Integrating over the super-cell 
volume we find that the total energy scales as ~ 1/ D^ ra 
where D max is the linear size of the super-cell in the z 
direction, i.e as ~ 1/n 3 . In Fig. (J^b) we report the BE 
obtained within the LRDMC method as a function of 
~ 1/n 3 . The BE curves close to the minimum for the 
32 and 64 atoms cells are reported in Fig. Extrap- 
olating the results reported in Fig. ([2b) to an infinite 
number of layers, we obtain a value of the BE of 60(5) 
meV/atom. This is reduced to 56(5), after adding zero 
point motion (~ 2 meV/atom) and lattice vibrational 
contributions at 300 K (~ 2 meV/atom), as computed 
from vibrational free energies using the data of Ref. [28| 
for phonon frequencies. 

Absorption experiments of aromatic molecules on 
graphite provide a measurement of the EE, while 
the cleavage energy is estimated to be 18% larger than 
exfoliation, on the basis of force field calculations [l2j ]. 
The cleavage energy is close, althogh not identical ot the 
BE defined here [29( . Therefore our comparison with ex- 
periment can only be indirect, as we computed BE, while 
experiment reports EE. Nevertheless it appears that our 
computed BE (~ 56 meV/atom ) is in good agreement 
with the estimate for the cleavage energy from the most 
recent experiments: 62 meV/atom. We note that in the 
analysis of experimental results, one makes use of fitted 
force fields to evaluate the contribution to EE of carbon- 
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mental measurements. The interaction energy between 
planes varies as D~ A 2 , i.e with a power law close to that 
found for two semiconducting planes. 
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FIG. 3: Binding Energy (BE) of graphite as obtained at the 
VMC (black solid circles) and LRDMC (red squares) level 
of theory, using a 2 x 2 x 1 super-cell, and at the LRDMC 
level (blue triangles) with a 2 x 2 x 2 super-cell. All data 
include finite size corrections using the scheme of Ref [23| . 
Dotted and solid lines are obtained with a fit with the func- 
tion aexp(SD) + b/D 4 . The experimental inter-layer distance 
[27j is shown by the dashed line. The large difference in BE 
between the 2x2x1 and 2x2x2 super-cell calculations 
arise from large spurious interactions between periodic images 
of planes in the case of the 2x2x1 cell. These spurious in- 
teractions are still present in the case of the 2x2x2 cell, 
but they are greatly reduced, as shown by the difference (~ 5 
meV/ atom) of BE obtained with 4 planes in the cell and the 
extrapolated value. 



hydrogen bonds. This is needed because an extrapolation 
is made on hydrocarbon adsorption energies, as a func- 
tion of the number of C atoms. Although the force field 
parameters were adjusted to experimental data, it is un- 
clear whether the assumption of additivity of forces close 
to the minimum is fully justified. In addition, the ratio 
between EE and BE is at present unknown and could 
only be estimated. 

While the value of the BE can be extrapolated for an 
infinite number of planes (as in Fig. ((2|)), the value of 
the equilibrium separation cannot. From the minima of 
the curves reported in Fig. ^ we obtain 3.350(24) and 
3.243(26) A at the VMC and LRDMC level, respectively, 
for the 2 x 2 x 1 cell, and 3.426(36) A at the LRDMC 
level, for the 2x2x2 cell. Difficulties arising from very 
fiat BE curves and, most importantly, from the lack of 
an extrapolation procedure as a function of the number 
of layers prevent us to find a fully converged equilibrium 
bond-length. The value found for the 2x2x2 cell is in 
good agreement with experiment (2% overestimate) (32|. 

In conclusion, we have investigated the bonding prop- 
erties of graphite in AB stacking, providing for the first 
time an estimate of the binding energy and long range be- 
havior of the total energy based on ab-initio, many body 
theory. Our calculated binding energy is in good agree- 
ment with most recent experiments, providing a bench- 
mark result for future calculations and further experi- 
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[29] The BE is defined as E n -i ay er — n * E\ — layer ; where 
E n -i a y er is the total energy of a graphite periodic ceil 
containing n layers at the optimized, equilibrium posi- 
tion and Ei-iayer is the total energy of a graphene layer 

[30] Corrections to the two-body term may also be treated 
using the method proposed in Ref. [25||, which requires 
the knowledge of the structure factor S(q) for q — » 
vectors. In principle 5*(q) can be directly evaluated using 
QMC calculations. In practice, for semi-metallic graphite 
and grap hene small error bars are difficult to achieve. In 
Ref. |26|1 the FS correction schemes proposed in Ref. [25l ] 
and [23| were shown to perform equally well. 



[31] We note that sums of pair- wise forces proportional to 
1/R 6 yield a behavior of the type E ~ D~ 4 . Of course 
our results showing a ~ D~ 4 behavior are not to be taken 
as confirming that description of binding in graphite with 
pair- wise 1/R e forces is correct. 

[32] For the 2x2x2 super-cell the position of the minimum 
was estimated considering LRDMC calculations at lattice 
space of a — 0.4 A . In the case of the 2x2x1 cell we 
observed that the extrapolation procedure for a — ► does 
not affect the equilibrium position. 



